SGD и оптимизация

SGD

Градиентный спуск: \[\theta_{k+1}=\theta_k-\eta\nabla f\left(\theta_{k}\right)\]

Стохастический градиентный спуск: \[\theta_{k+1}=\theta_k-\eta\hat{\nabla f\left(\theta_{k}\right)}\]

Осуществление шага

Рассмотрим несколько методик осуществления шага и получения ответа

• Обычный стохастический градиентный спуск: \[\theta_{k+1}=\theta_k-\eta\hat{\nabla f\left(\theta_{k}\right)}\]
• Nesterov Accelerated Gradient Descent: \[ y_{k+1}=\theta_k-\eta\hat{\nabla f\left(\theta_k\right)}\] \[ \theta_{k+1}=\left(1-\gamma\right)y_{k+1}+\gamma y_k\]
• SGD с моментом: \[ v_0 = \mathbf{0}\] \[ v_{k+1} = \alpha v_{k} - \eta\hat{\nabla f\left(\theta_k\right)}\] \[ \theta_{k+1} = \theta_k+v\]
• SGD с усреднением: \[ \hat{\theta} = \frac{1}{n}\sum\limits_{k=1}^{n}\theta_k\]
• Неявный SGD: \[ \theta_{k+1}=\theta_k-\eta\hat{\nabla f\left(\theta_{k+1}\right)}\]

Обновление learning rate:

• \(\eta_k=\eta_0\left(1+\alpha k\right)^{-c}\)
• AdaGrad: \[ G_k=\sum\limits_{i=1}^k \hat{\nabla f\left(\theta_k\right)}\hat{\nabla f\left(\theta_k\right)}^\mathrm{T} \] \[ \theta_{k+1}=\theta_k-\eta \mathrm{diag}\left(G_k\right)^{-\frac{1}{2}} \hat{\nabla f\left(\theta_k\right)}\]
• RMSProp: \[ v_{k+1}=\gamma v_k+\left(1-\gamma\right)\mathrm{diag}\left( \hat{\nabla f\left(\theta_k\right)}\hat{\nabla f\left(\theta_k\right)}^\mathrm{T}\right)\] \[ \theta_{k+1}=\theta_k-\eta \mathrm{diag}\left(v_{k+1}\right)^{-\frac{1}{2}}\hat{\nabla f\left(\theta_k\right)}\]

SGD в R

Для демонстрации SGD в R будем использовать функцию sgd из пакета sgd

sgd(formula, data, model, model.control = list(), sgd.control = list(...))

• formula – описание зависимости
• data – данные
• model – используемая модель: “lm”, “glm”, “cox” (Cox proportional hazards model), “gmm” (generalized method of moments), “m” (m-estimation)
• model.control – параметры модели
• sgd.control – список параметров процедуры SGD:
  • method – используемый метод SGD: “sgd”, “implicit”, “asgd”, “ai-sgd”, “momentum”, “nesterov”
  • lr – процедура обновления learning rate: “one-dim”, “one-dim-eigen”, “d-dim”, “adagrad”, “rmsprop”
  • lr.control – параметры выбраной процедуры обновления learning rate
  • start – начальный вектор параметров
  • reltol – критерий остановки (в смысле относительного изменения параметров)
  • npasses – количество проходов по данным
  • shuffle – перемешивание данных

Возвращаемый объект будет содержать поля

• converged – индикатор сходимости
• pos – моменты времени, в которые сохранено внутреннее состояние
• estimates – оценки параметров в моменты pos
• times – время, затраченное на каждую из итераций
• coefficients – оценки коэффициентов

Функция позволяет оценивать параметры следующих моделей:

• Линейная и обобщённая линейная модели:
  • formula – формула описывающая линейную модель
  • model.control – список:
    • family – семейство (как в функциии glm; по-умолчанию – gaussian)
• Обобщённый метод моментов:
  • x – выборка
  • model.control – список:
    • fn – функция, возращающая условия моментов
    • gr – функция, возращающая градиенты условий моментов (иначе они будут вычисляться из fn числено)
    • nparams – число параметров
• M-оценки:
  • formula – формула, описывающая линейную модель
  • model.control – список:
    • loss – функция потерь (на текущий момент поддерживается только функция потерь Huber’а)
• Модель Cox’а:
  • formula – формула, описывающая модель пропорциональных рисков

Для сравнения расширений SGD построим линейную модель \(y=x_1\beta_1+x_2\beta_2+\varepsilon\) и рассмотрим изменение оценок параметров во времени:

set.seed(42)
N <- 10000
sigma <- 0.1
beta <- c(5, 1)
start <- c(-15, -15)
x.1 <- rnorm(N)
x.2 <- rnorm(N)
y <- beta[1] * x.1 + beta[2] * x.2 + rnorm(N) * sigma
df <- data.frame(x.1=x.1, x.2=x.2, y=y)

methods <- c('sgd', 'nesterov', 'ai-sgd', 'implicit', 'momentum')
models <- lapply(methods, function(m)sgd(y~x.1+x.2+0, data=df, model='lm', sgd.control=list(method=m, start=start, shuffle=F, npasses=3, pass=T)))
dd <- do.call(rbind, lapply(1:length(models), function(mi){m<-models[[mi]]; data.frame(xy=rbind(start, t(m$estimates)), method=rep(methods[mi], length(m$pos)+1))}))
colnames(dd) <- c('beta_1', 'beta_2', 'method')

Пути, полученные различными алгоритмами и линии уровня целевой функции

x.min <- -20
x.max <- 20
xy <- meshgrid(seq(x.min, x.max, by=0.1))
X <- as.vector(xy$X)
Y <- as.vector(xy$Y)
Z <- sapply(1:length(X), function(id)sum((X[id]*x.1+Y[id]*x.2-y)^2/length(x.1)))
h<-xyplot(beta_1~beta_2, data=dd, type='l', groups=dd$method, auto.key=T, scales=list(x=list(limits=c(x.min, x.max)), y=list(limits=c(x.min, x.max))), main='SGD paths & target function levels')
h<-h+contourplot(Z~Y*X, region=F,col.regions=brewer.pal(11,'Spectral'), colorkey=F)
h<-h+xyplot(beta[1]~beta[2], pch=4, cex=2, col='black')
h<-h+xyplot(start[1]~start[2], pch=3, cex=2, col='black')
h

coeffs<-lapply(1:length(methods), function(id)list(method=methods[id], coeffs = models[[id]]$coefficients))

coeffs

## [[1]]
## [[1]]$method
## [1] "sgd"
## 
## [[1]]$coeffs
##          [,1]
## [1,] 5.000435
## [2,] 1.001512
## 
## 
## [[2]]
## [[2]]$method
## [1] "nesterov"
## 
## [[2]]$coeffs
##          [,1]
## [1,] 4.999665
## [2,] 1.001213
## 
## 
## [[3]]
## [[3]]$method
## [1] "ai-sgd"
## 
## [[3]]$coeffs
##          [,1]
## [1,] 4.998543
## [2,] 1.000185
## 
## 
## [[4]]
## [[4]]$method
## [1] "implicit"
## 
## [[4]]$coeffs
##          [,1]
## [1,] 4.999301
## [2,] 1.000442
## 
## 
## [[5]]
## [[5]]$method
## [1] "momentum"
## 
## [[5]]$coeffs
##          [,1]
## [1,] 4.999665
## [2,] 1.001213

Для оптимизации произвольных функций методом стохастического градиентного спуска можно воспользоваться обобщённым методом моментов, определив градиент функции;

Например, для функции Розенброка (\(f\left(x, y\right)=\left(1-x\right)^2+100\left(y-x^2\right)^2\)):

set.seed(42)
N <- 10000
start <- c(-1.5, 3)
beta <- c(1, 1)
sigma <- 0.1
x<-matrix(rnorm(N*2), ncol=2)*sigma
y<-as.matrix(rep(NA, N), ncol=1)
gr<-function(theta, xx) {
    x<-theta[1]
    y<-theta[2]
    return(as.matrix(c(400*x^3-400*x*y+2*x-2+xx[1], 200*(y-x^2)+xx[2])))
}

methods <- c('sgd')
models <- lapply(methods, function(m)
    sgd(x, y=y, model='gmm', model.control=list(gr=gr, nparams=2), sgd.control=list(method=m, start=start, shuffle=F, npasses=6, pass=T, lr='adagrad')))

dd <- do.call(rbind, lapply(1:length(models), function(mi){m<-models[[mi]]; data.frame(xy=rbind(start, t(m$estimates)), method=rep(methods[mi], length(m$pos)+1))}))
colnames(dd) <- c('beta_1', 'beta_2', 'method')
x.min <- -2
x.max <- 4
xy <- meshgrid(seq(x.min, x.max, by=0.1))
X <- as.vector(xy$X)
Y <- as.vector(xy$Y)
Z <- sapply(1:length(X), function(id){x=X[id]; y=Y[id]; (1-x)^2+100*(y-x^2)^2})
h<-xyplot(beta_1~beta_2, data=dd, type='l', groups=dd$method, auto.key=T, scales=list(x=list(limits=c(x.min, x.max)), y=list(limits=c(x.min, x.max))), main='SGD paths & target function levels')
h<-h+contourplot(Z~Y*X, region=F,at=(1:10)*10,col.regions=brewer.pal(11,'Spectral'), colorkey=F)
h<-h+xyplot(beta[1]~beta[2], pch=4, cex=2, col='black')
h<-h+xyplot(start[1]~start[2], pch=3, cex=2, col='black')
h

coeffs<-lapply(1:length(methods), function(id)list(method=methods[id], coeffs = models[[id]]$coefficients))

coeffs

## [[1]]
## [[1]]$method
## [1] "sgd"
## 
## [[1]]$coeffs
##           [,1]
## [1,] 0.9997515
## [2,] 0.9992494

Зависимость от learning rate

Сгенерируем новый набор данных с большим количеством параметров и рассмотрим зависимость ошибки от итерации при различном выборе начального learning rate

Сравним методику с адаптивным изменением learning rate и степенным убыванием

set.seed(42)
N <- 10000
p <- 40
sigma <- 0.1
tform <- diag(runif(p, min=-1, max=1)*20) #%*% matrix(runif(p*p), nrow=p)
data <- matrix(rnorm(p * N), nrow = N)
data <- data %*% tform

beta <- rnorm(p)
beta.0 <- rnorm(1)
reference <- c(beta.0, beta)
y <- data %*% beta + beta.0  + rnorm(N) * sigma

data <- data.frame(y=y, x=data)

rates <- 10^((-2:4))

sgd.models <- lapply(rates, function(r)sgd(y~., data=data, model='lm', sgd.control=list( lr.control=c(r, NA, NA, NA), pass=T, npasses=50, start=rep(0, length(beta)+1))))
sgd.models.rmsprop <- lapply(rates, function(r)sgd(y~., data=data, model='lm', sgd.control=list(lr='rmsprop', lr.control=c(r, NA, NA), pass=T, npasses=50, start=rep(0, length(beta)+1))))

Построим зависимость конечной ошибки в оцениваемых параметрах от learning rate

errors <- sapply(sgd.models, function(m)norm(m$coefficients - reference))
errors.rmsprop <- sapply(sgd.models.rmsprop, function(m)norm(m$coefficients - reference))
df <- data.frame(rates=rates, errors.sgd=errors, errors.rmsprop=errors.rmsprop)
df

##   rates errors.sgd errors.rmsprop
## 1 1e-02 1.09077874     0.07253628
## 2 1e-01 0.11802718     0.06313638
## 3 1e+00 0.06529551     0.06214990
## 4 1e+01 0.06230022     0.06205905
## 5 1e+02 0.06207420     0.06204986
## 6 1e+03 0.06205138     0.06204894
## 7 1e+04 0.06204910     0.06204885

Можно сделать вывод, что использование специальной методики обновления learning rate позволяет практически исключить зависимость результата от правильного выбора его начального значения.

Ниже приведены зависимость ошибки оценки на очередном шаге для степенного убывания learning rate:

lines <- lapply(sgd.models,  function(m)data.frame(pos=m$pos, error=apply(m$estimates, 2, FUN=function(r)norm(as.matrix(r - reference)))))

lines.lr <- lapply(1:length(rates), function(i){l <- lines[[i]]; lr <- rep(rates[i], dim(l)[1]); data.frame(l, learning.rate=lr)})
lines.all <- do.call(rbind, lines.lr)
xyplot(error~pos, lines.all, groups=lines.all$learning.rate, type='l', auto.key=T, scales=list(y=list(log=T)), main='Error vs iteration number: decreasing learning rate')

Ниже приведены зависимость ошибки оценки на очередном шаге для обновления learning rate с помощью rmsprop:

lines.rmsprop <- lapply(sgd.models.rmsprop,  function(m)data.frame(pos=m$pos, error=apply(m$estimates, 2, FUN=function(r)norm(as.matrix(r - reference)))))

lines.rmsprop.lr <- lapply(1:length(rates), function(i){l <- lines.rmsprop[[i]]; lr <- rep(rates[i], dim(l)[1]); data.frame(l, learning.rate=lr)})
lines.rmsprop.all <- do.call(rbind, lines.rmsprop.lr)
xyplot(error~pos, lines.rmsprop.all, groups=lines.rmsprop.all$learning.rate, type='l', auto.key=T, scales=list(y=list(log=T)), main='Error vs iteration number: adaptive learning rate')

//: //: Слишком быстрое убывание learning rate может //: {r,warning=F,fig.width=10,fig.height=10} [//]: set.seed(42) [//]: sigma <- 0.1 [//]: x <- seq(-1, 1, by=0.01) [//]: y <- c(abs(x)) [//]: data <- data.frame(x=x, y=y+rnorm(length(x))*sigma) [//]: [//]: rates <- 10^((-5:5)) [//]: [//]: sgd.models <- lapply(rates, function(r)sgd(y~poly(x, 2), data=data, model='lm', sgd.control=list( lr.control=c(r, NA, NA), pass=T, npasses=100, lr='rmsprop'))) [//]: [//]: lm(y~poly(x, 2)) [//]: sgd.models [//]:

Скорость сходимости в сравнении с OLS

Сравним временную сложность SGD с OLS

Так как для задачи с \(N\) наблюдениями и \(p\) параметрами вычислительная сложность имеет порядок \(O\left(Np^2\right)\), для демонстрации превосходства SGD выберем количество параметров достаточно большим.

set.seed(42)
N <- 20000
p <- 1000
sigma <- 1
data <- matrix(rnorm(p * N), nrow = N)

beta <- rep(5, p)
beta.0 <- 5
reference <- c(beta.0, beta)
y <- data %*% beta + beta.0  + rnorm(N) * sigma

data <- data.frame(y=y, x=data)

rates <- 1

time.sgd<-system.time(sgd.models.rmsprop <- sgd(y~., data=data, model='lm',  npasses=3, reltol=1e-9, lr='rmsprop', start=rep(0, length(beta)+1), size=100))
time.lm<-system.time(lm.model <- lm(y~., data=data))

residuals.sgd <- as.matrix(data[,-1])%*%sgd.models.rmsprop$coefficients[-1]+sgd.models.rmsprop$coefficients[1]-y

m <- sgd.models.rmsprop
errors.rmsprop <- c(norm(residuals.sgd)/norm(y), time.sgd[3])
m <- lm.model
errors.lm <- c(norm(as.matrix(m$residuals))/norm(y), time.lm[3])
errors <- rbind(errors.rmsprop, errors.lm)
colnames(errors) <- c('error',  'time')
rownames(errors) <- c('SGD',  'OLS')

errors

##          error   time
## SGD 0.02387713  6.029
## OLS 0.00619758 21.381

Заметим также, что для задачи нелинейной регрессии время, требуемое в SGD на вычисление величины шага не меняется, в то время как нелинейная регрессия перестаёт иметь решение в явном виде (и требует применения итеративной процедуры).

LARS и lasso с помощью LARS

Продемонстриреум отбор признаков с помощью LARS на примере набора данных diabetes (используя пакет lars)

Датасет содержит в себе 10 оригинальных признаков и набор взаимодействий между ними.

data <- as.data.frame(diabetes)
corrplot(cor(data))

Совершим регрессию с помощью LARS:

lars(x, y, type = c("lasso", "lar", "forward.stagewise", "stepwise", trace = FALSE, normalize = TRUE, intercept = TRUE, Gram, eps = .Machine$double.eps, max.steps, use.Gram = TRUE)

• x – признаки
• y – предсказываемая величина
• type – метод (lasso сводится к lar с ограничениями)
• trace – вывод прогресса
• normalize – стандартизация входных значений
• intercept – добавление свободного члена
• Gram – матрица корреляций (для повторных запусков)
• eps – машинный \(\varepsilon\)
• max.steps – число шагов

Для того, чтобы выбрать число признаков, воспользуемся кросс-валидацией

cv.lars(x, y, K = 10, index, trace = FALSE, plot.it = TRUE, se = TRUE, type = c("lasso", "lar", "forward.stagewise", "stepwise"), mode = c("fraction", "step"), ...)

• x – признаки
• y – предсказываемая величина
• type – метод (lasso сводится к lar с ограничениями)
• K – кратность кросс-валидации
• index – момент (регуляризации), в который вычисляются значения коэффициентов и оценка ошибки по кросс-валидации
• trace – вывод прогресса
• _plot.it__ – построение графика
• se – строить оценку стандартной ошибки
• mode – режим кросс-валидации:
  • fraction – доля от насыщенной модели (для lasso и forward.stagewise)
  • step – по шагам (для stepwise и lar)

X <- as.matrix(data[,colnames(data)!='y'])
cv.result <- cv.lars(X, data$y, index=1:50, type='lar', mode='step')

Можно видеть, что минимум достигается на количестве элементов

min.step = cv.result$index[which.min(cv.result$cv)]
print(min.step)

## [1] 16

Оценим коэффициенты модели

model.lars <- lars(X, data$y, type='lar', max.steps=50)
print(model.lars$beta[min.step, model.lars$beta[min.step,]!=0])

##       x.sex       x.bmi       x.map       x.glu      x2.sex      x2.hdl 
##  -94.904160  501.616811  251.792232   18.130314  -17.408382 -187.827872 
##      x2.ltg    x2.age^2    x2.bmi^2    x2.glu^2  x2.age:sex  x2.age:map 
##  467.778725    7.608734   38.748381   69.809696  107.727556   30.045318 
##  x2.age:ltg  x2.age:glu  x2.bmi:map 
##    8.613394   11.600284   85.688057

Также можно построить изменение оценок коэффициентов по шагам:

plot(model.lars)